#setwd('/data1/bsi/bioinf_int/s106381.borawork/data_HTHGU/HTHGU_GeneSumm/samplefiles_out')
args <- commandArgs(TRUE)
setwd(args[1])
#library(genefilter, lib.loc=args[2])
#library(gcrma, lib.loc=args[2])
#library(Biobase, lib.loc=args[2])
#library(preprocessCore, lib.loc=args[2])
#library(MASS, lib.loc=args[2])
#library(affy, lib.loc=args[2])
#library(simpleaffy, lib.loc=args[2])
#library(frma, lib.loc=args[2])
#library(hgu133plus2geneccdscdf, lib.loc=args[2])
#library(hgu133plus2frmavecs,lib.loc=args[2])
#library(hgu133afrmavecs, lib.loc=args[2])
#library(equalprobesorihgu133ageneccdscdf, lib.loc=args[2])
#library(affyPLM, lib.loc=args[2])


.libPaths("/usr/local/biotools/r/R-2.14.0/lib64/R/library")
library(genefilter)
library(gcrma)
library(Biobase)
library(preprocessCore)
library(MASS)
library(affy)
library(AnnotationDbi)
library(affyPLM)
#library(simpleaffy, lib=args[2])
.libPaths("/home/m081429/R/x86_64-unknown-linux-gnu-library/2.12")
library(frma)
library(hgu133plus2frmavecs)
#library(hgu133ahsentrezgcdf, lib=args[2])
library(hgu133plus2hsentrezgcdf)


data <- ReadAffy()
#data@cdfName<-"hgu133plus2geneccds"
data@cdfName<-"hgu133plus2hsentrezg"
#data@cdfName<-"equalprobesorihgu133ageneccds"
#script start
#######################################################################################
frma1 <- function(object, background="rma", normalize="quantile", summarize="robust_weighted_average", target="core", input.vecs=list(normVec=NULL, probeVec=NULL, probeVarBetween=NULL, probeVarWithin=NULL, probesetSD=NULL), output.param=NULL, verbose=FALSE){

  if(!class(object) %in% c("AffyBatch", "ExonFeatureSet")) stop("object must be of class AffyBatch or ExonFeatureSet.")

  if(class(object)=="ExonFeatureSet") target <- match.arg(target, c("core", "full", "extended", "probeset"))
  if(verbose & class(object)=="ExonFeatureSet") message(paste("Exon summarization at ", target, " level.\n", sep=""))

  if(!background %in% c("none", "rma")) stop("background must be either none or rma")
  if(!normalize %in% c("none", "quantile")) stop("normalize must be either none or quantile")
  if(!summarize %in% c("median_polish", "average", "median", "weighted_average", "robust_weighted_average", "random_effect")) stop("summarize must be one of: median_polish, average, median, weighted_average, robust_weighted_average, random_effect")

  if(summarize=="random_effect" & (any(input.vecs$probeVarBetween==0) | any(input.vecs$probeVarWithin==0))) stop("If summarize method is random_effect then probeVarBetween and probeVarWithin must be greater than zero for all probes.")

  if(summarize=="random_effect" & class(object)=="ExonFeatureSet") stop("Summarization method random_effect is not implemented for ExonFeatureSet objects.")
  
  if(class(object)=="AffyBatch"){
    cdfname <- cleancdfname(cdfName(object))
    platform <- gsub("cdf","",cdfname)
  }
  if(class(object)=="ExonFeatureSet"){
    cdfname <- annotation(object)
    platform <- paste(cdfname, target, sep="")
  }

  if(summarize == "median_polish") output <- frmaMedPol(object, background, normalize, target, input.vecs, verbose)
  if(summarize %in% c("average", "median", "weighted_average", "robust_weighted_average")) output <- frmaRobReg(object, background, normalize, summarize, target, input.vecs, output.param, verbose)
  if(summarize == "random_effect") output <- frmaBatch1(object, background, normalize, input.vecs, output.param, verbose)

  if(is.null(output.param)){
    if(is.null(output$stderr)){
      e <- new("ExpressionSet", assayData=assayDataNew(exprs=output$exprs), annotation=platform)
    } else{
      e <- new("ExpressionSet", assayData=assayDataNew(exprs=output$exprs, se.exprs=output$stderr), annotation=platform)
    }
  } else {
    if(is.null(output$stderr)){
      e <- new("frmaExpressionSet", assayData=assayDataNew(exprs=output$exprs), annotation=platform)
    } else{
      e <- new("frmaExpressionSet", assayData=assayDataNew(exprs=output$exprs, se.exprs=output$stderr), annotation=platform)
    }
    if("weights" %in% output.param) w <- output$weights else w <- matrix(nrow=0, ncol=0)
    if("residuals" %in% output.param) r <- output$residuals else r <- matrix(nrow=0, ncol=0)
    if("random_effects" %in% output.param) g <- output$gammas else g <- matrix(nrow=0, ncol=0)
    e@weights <- w
    e@residuals <- r
    e@randomeffects <- g
  }

  return(e)
}

########################################################################################


#########################################################################################
frmaBatch1 <- function(object, background, normalize, input.vecs, output.param, verbose){

cdfname <- cleancdfname(cdfName(object))
cdfname<-"hgu133plus2"


  tmp <- gsub("cdf","",cdfname)
  platform <- gsub("..entrezg", "", tmp)
  vecdataname <- paste(tmp, "frmavecs", sep="")
  
  if(background == "rma"){
    if(verbose) message("Background Correcting ...\n")
    object <- bg.correct.rma(object)
  }

  if(is.null(input.vecs$normVec) | is.null(input.vecs$probeVec) | is.null(input.vecs$probeVarWithin) | is.null(input.vecs$probeVarBetween)){
    pkg <- paste(platform, "frmavecs", sep="")
    require(pkg, character.only=TRUE, quietly=TRUE) || stop(paste(pkg, "package must be installed first"))
    data(list=eval(vecdataname))

    if(is.null(input.vecs$normVec)) input.vecs$normVec <- get(vecdataname)$normVec
    if(is.null(input.vecs$probeVec)) input.vecs$probeVec <- get(vecdataname)$probeVec
    if(is.null(input.vecs$probeVarWithin)) input.vecs$probeVarWithin <- get(vecdataname)$probeVarWithin
    if(is.null(input.vecs$probeVarBetween)) input.vecs$probeVarBetween <- get(vecdataname)$probeVarBetween
  }

  if(normalize == "quantile"){
    if(verbose) message("Normalizing ...\n")
    pm(object) <- normalize.quantiles.use.target(pm(object), input.vecs$normVec)
  }

  if(verbose) message("Summarizing ...\n")
  pns <- probeNames(object)
  pms <- log2(pm(object))

  N <- 1:dim(pms)[1]
  S <- split(N, pns)

  fit <- lapply(1:length(S), function(i) {
    s <- S[[i]]
    batchFit(pms[s,], input.vecs$probeVec[s], input.vecs$probeVarWithin[s], input.vecs$probeVarBetween[s])
  })
  
  exprs <- matrix(unlist(lapply(fit, function(x) x$exprs)), ncol=ncol(pms), byrow=TRUE)
  rownames(exprs) <- names(fit)
  colnames(exprs) <- colnames(pms)

  if("weights" %in% output.param){
    weights <- matrix(unlist(lapply(fit, function(x) x$w)), ncol=ncol(pms), byrow=TRUE)
    rownames(weights) <- pns
    colnames(weights) <- sampleNames(object)
  } else weights <- NULL

  if("residuals" %in% output.param){
    residuals <- apply(exprs,2,function(x) rep(x, table(pns)))
    residuals <- (pms-input.vecs$probeVec) - residuals
    rownames(residuals) <- pns
    colnames(residuals) <- sampleNames(object)
  } else residuals <- NULL

  if("random_effects" %in% output.param){
    gammas <- matrix(unlist(lapply(fit, function(x) x$gamma)), ncol=1)
    rownames(gammas) <- pns
    colnames(gammas) <- NULL
  } else gammas <- NULL
  
  colnames(exprs) <- sampleNames(object)
  rownames(exprs) <- unique(pns)
  
  return(list(exprs=exprs, stderr=NULL, weights=weights, residuals=residuals, gammas=gammas))
}

batchFit <- function(y, probe, sigma, tau, maxiter=100){
  nchip <- ncol(y)
  nprobe <- nrow(y)
  y.new <- as.vector(y - probe)  #kill any names
  n <- length(y)
  wt <- matrix(rep(1/as.vector(sigma), nchip), nrow=nprobe)

  #rlm setup, taken from the rlm function
  irls.delta <- function(old, new)
      sqrt(sum((old - new)^2)/max(1e-20, sum(old^2)))
  psi.huber <- function(u, k=1.145)  pmin(1, k/abs(u))

  # Setup for linear regression
  xtx <- matrix(0., ncol=(nchip+nprobe), nrow=nchip+nprobe)
  i1 <- seq(1, length=nchip+nprobe, by=nchip+nprobe+1) #index of diagonal
  i2a <- which(row(xtx) <= nprobe & col(xtx) >nprobe)
  i2b <- c(outer(seq(nprobe, length=nprobe, by=nrow(xtx)), 1:nchip,  '+'))
  i3 <- as.vector(row(y))          #gamma coefs
  i4 <- as.vector(col(y)) + nprobe #theta coefs
  sqwt <- sqrt(wt) 
  w <- rep(1., length(y))

  # Initial linear regression
  xtx[i1] <- c(rowSums(wt) + 1/tau, colSums(wt))
  xtx[i2a] <- wt
  xtx[i2b] <- wt         
  xty <- c(rowSums(wt*y.new), colSums(wt*y.new))

  beta <- solve(xtx, xty)
  resid <- sqwt * (beta[i3] + beta[i4] - y.new)  # weighted residuals

  for (iter in 1:maxiter) {
    oldresid <- resid
    scale <- median(abs(resid))/.6745
    w <- matrix(psi.huber(resid/scale, k=1.345), nrow=nprobe)
    w2 <- w*wt
    xtx[i1] <- c(rowSums(w2) + 1/tau, colSums(w2))
    xtx[i2a] <- w2
    xtx[i2b] <- w2
    temp <- w2* y.new
    xty <- c(rowSums(temp), colSums(temp))
    
    beta <- solve(xtx, xty)
    resid <- sqwt * (beta[i3] + beta[i4] - y.new)  # weighted residuals
    convi <- irls.delta(oldresid, resid)
    if (convi < 1e-4) break
  }
  return(list(exprs= beta[-(1:nprobe)], w=as.vector(t(w)), gamma=beta[1:nprobe]))
}



#script end

affy.dc5 <- mas5calls(data)
affy.pval <- assayData(affy.dc5)[["se.exprs"]]
#bc <- barcode(data)
#pval <- detection.p.val(data)
#pval.c <- pval$call
#pval.p <- pval$pval
#pval.pc <- cbind(pval.p,pval.c)
object <- frma1(data, summarize = "random_effect")
mydata <- as.vector(strsplit(args[1], '/')[[1]])
output.file=paste(mydata[length(mydata)],".normalized",sep="")
output.file.call = paste(output.file,".call",sep="")
output.file.barcode = paste(output.file,".barcode",sep="")
write.exprs(object,output.file,sep="\t")
#write.table(coefs(object),output.file,sep="\t",quote=FALSE)
write.table(affy.pval,output.file.call,sep="\t",quote=FALSE)
#write.table(bc,output.file.barcode,sep="\t",quote=FALSE)
#dat<-read.table(output.file,sep="\t",header=T)
#final_dat<-cbind(dat,pval.pc)
#nam <-as.vector(colnames(final_dat))
#nam[1]<-"Probeset"
#nam[5]<-"Sample Name"
#nam[6]<-"Affy_Platform"
#final_dat<-cbind(final_dat,nam[2],annotation(data))
#colnames(final_dat)<-nam
#write.table(final_dat,output.file,sep="\t",row.names=F,quote=F)
q()
n
